**************************************************
* Author: Michael Tomz
* Date: Dec 20, 2001
* Format: Stata do-file
* Purpose: Replicate Section 5.1 (Normal vs. T)
**************************************************

version 6.0
clear
set mem 40m

* LOAD THE RESULTS
use r92, clear
append using r87
append using r83
append using r79
append using r74a
append using r74
append using r70
append using r66
append using r64
append using r59


* ARE SUR ERRORS (FOR FULLY CONTESTED DISTRICTS) NORMALLY DISTRIBUTED?
* Create a matrix called `cover' that shows what percent of residuals 
* fell within the 95%, 90%, 80%, and 50% elliptical confidence regions.
* Did 95% of residuals fall within the 95% region, as expected?  Etc.
tempname cover                               /* ci coverage */
matrix `cover' = J(11,4,0)                   /* matrix to hold coverage values */
matrix colnames `cover' = 95pct 90pct 80pct 50pct
matrix rownames `cover' = 59 64 66 70 74 74a 79 83 87 92 avg
local i 1
local years 59 64 66 70 74 74a 79 83 87 92
while "`years'" ~= "" {
   quietly {
      gettoken year years : years            /* loop through the years */
      preserve
      keep if yearst == "`year'" & fully==1  /* analyze fully contested districts for each year */
      gen err1c = y1 - ey1c                  /* error using clarify's SUR procedure, equation 1 */
      gen err2c = y2 - ey2c                  /* error using clarify's SUR procedure, equation 2 */
      count if err1c~=. & err2c~=. & yearst == "`year'"
      local numobs = r(N)
      sureg (y1 y1lag y2lag coninc labinc libinc) (y2 y1lag y2lag coninc labinc libinc) if yearst == "`year'"
      mat Sig = e(Sigma)
      mat Siginv = syminv(Sig)
      * Calculate Mahalanobis distance
      gen dist = err1c*(err1c*Siginv[1,1]+err2c*Siginv[1,2]) + err2c*(err1c*Siginv[1,2]+err2c*Siginv[2,2])
      * Compare to quantiles of Chi-Squared distribution with 2 d.f.
      count if dist <= invchi(2,.05)  /* 95% contour */
      matrix `cover'[`i',1] = 100*r(N)/`numobs'
      matrix `cover'[11,1] = `cover'[11,1] + 10*r(N)/`numobs'
      count if dist <= invchi(2,.1)   /* 90% contour */
      matrix `cover'[`i',2] = 100*r(N)/`numobs'
      matrix `cover'[11,2] = `cover'[11,2] + 10*r(N)/`numobs'
      count if dist <= invchi(2,.2)   /* 80% contour */
      matrix `cover'[`i',3] = 100*r(N)/`numobs'
      matrix `cover'[11,3] = `cover'[11,3] + 10*r(N)/`numobs'
      count if dist <= invchi(2,.5)   /* 50% contour */
      matrix `cover'[`i',4] = 100*r(N)/`numobs'
      matrix `cover'[11,4] = `cover'[11,4] + 10*r(N)/`numobs'
      restore
   }
   local i = `i' + 1
}
di "Are SUR errors Normally Distributed?"
mat list `cover'

* ARE KK ERRORS T-DISTRIBUTED?
matrix `cover' = J(11,4,0)
matrix colnames `cover' = 95pct 90pct 80pct 50pct
matrix rownames `cover' = 59 64 66 70 74 74a 79 83 87 92 avg
local i 1
local years 59 64 66 70 74 74a 79 83 87 92
while "`years'" ~= "" {
   quietly {
      gettoken year years : years
      preserve
      keep if yearst == "`year'" & fully==1
      gen err1 = y1 - ey1
      gen err2 = y2 - ey2
      count if err1~=. & err2~=.
      local numobs = r(N)
      scalar df = tdf[1]                              /* pluck off the degrees of freedom */
      matrix Sig = (ts11[1],ts12[1]\ts12[1],ts22[1])  /* pluck off elements of Sigma for T distrib */
      * Under KK, Sigma is for the Normal, so no need to downscale with Sig = Sig*(df-2)/df
      mat Siginv = syminv(Sig)                       
      gen dist = err1*(err1*Siginv[1,1]+err2*Siginv[1,2]) + err2*(err1*Siginv[1,2]+err2*Siginv[2,2])
      count if dist <= invfprob(2,df,.05)*2   /* 95% Ci */
      matrix `cover'[`i',1] = 100*r(N)/`numobs'
      matrix `cover'[11,1] = `cover'[11,1] + 10*r(N)/`numobs'
      count if dist <= invfprob(2,df,.1)*2    /* 90% Ci */
      matrix `cover'[`i',2] = 100*r(N)/`numobs'
      matrix `cover'[11,2] = `cover'[11,2] + 10*r(N)/`numobs'
      count if dist <= invfprob(2,df,.2)*2    /* 80% Ci */
      matrix `cover'[`i',3] = 100*r(N)/`numobs'
      matrix `cover'[11,3] = `cover'[11,3] + 10*r(N)/`numobs'
      count if dist <= invfprob(2,df,.5)*2    /* 50% Ci */
      matrix `cover'[`i',4] = 100*r(N)/`numobs'
      matrix `cover'[11,4] = `cover'[11,4] + 10*r(N)/`numobs'
   }
   restore
   local i = `i' + 1
}
di "Are KK errors t-Distributed?"
mat list `cover'
